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SUMMARY 


The  following  sections  describe  in  report  form  for  the  first  time  the  HULL 
code.  Since  the  code  itself  is  quite  large  and  required  the  work  of  dozens 
of  people  to  produce  the  capability  that  the  code  now  offers,  one  cannot  expect 
to  learn  all  the  pertinent  details  from  this  report.  However,  one  can  begin  to 
appreciate  the  power  that  the  code  provides  as  a  scientific  tool.  The  basic 
HULL  code  is  a  hydrodynamics  code  that  numerically  solves  the  equations  of 
mass,  momentum  and  energy  conservation. 

How  these  equations  are  solved  and  how  the  results  compare  to  experimental 
evidence  will  be  examined. 
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PREFACE 

This  report  was  prepared  from  contributions  of  many  individuals  and  is 
meant  to  serve  as  an  introduction  to  the  HULL  computer  code.  The  development 
of  the  HULL  code  involved  a  large  number  of  people  and  involved  many  man-years 
of  effort.  We  would  like  to  express  our  thanks  to  the  following  people  who 
have  contributed  many  hours  of  hard  work  to  this  effort:  Reginald  Clemens  of 
Science  Applications  Inc  ,  Lewis  Gaby,  John  Prentice  and  Mike  Tower  all  from 
the  Air  Force  Weapons  Laboratory. 
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SECTION  I 
INTRODUCTION 

This  report  describes  the  HULL  code,  which  is  used  at  the  Air  Force  Weapons 
Laboratory  (AFWL)  to  theoretically  investigate  hydrodynamic  phenomena  of 
interest  to  the  Air  Force.  Usual  problems  of  interest  involve  atmospheric 
response  (e.g.,  air  blast,  heave,  etc.)  to  large-scale  rapid  energy  deposition 
(e.g.,  nuclear  and  high  explosive  detonations). 

The  HULL  code  (or  simply  HULL)  is  actually  a  system  of  computer  codes  used 
primarily  for  hydrodynamic  calculations,  and  hence  is  normally  referred  to  as 
a  hydrocode.  Where  possible,  its  operation  is  automatic  to  allow  the  user  time 
to  run  and  analyze  calculations  instead  of  doing  burdensome  bookkeeping.  The 
system  is  written  in  a  higher-level  language,  ..nich  when  compiled  generates 
Control  Data  Corporation  (COC)  extended  FORTRAN  and  a  few  somewhat  system- 
oriented  routines  written  in  assembly  language.  Utilizing  the  higher-level 
language  allows  one  to  produce  a  large  variety  of  codes  by  specifying  certain 
options  at  the  beginning  of  a  run. 

The  HULL  code,  initiated  in  September  1971,  had  evolved  through  88  separate 
versions  by  June  1976.  It  continues  to  be  modified  in  order  to  allow  the 
addition  of  new  options  as  well  as  improvements  in  numerical  techniques.  The 
current  version  consists  of  approximately  44000  lines  of  code  and  can  be 
obtained  from  the  AFWL  for  users  with  legitimate  needs  for  the  code. 

Basically,  HULL  is  a  multidimensional  (two  or  three)  Eulerian  hydrodynamics 
code.  For  two-dimensional  problems,  Cartesian  or  cylindrical  geometry  is  avail¬ 
able,  whereas  for  three-dimensional  problems  only  Cartesian  geometry  is  avail¬ 
able  at  present.  Two  finite  difference  schemes  may  be  used— the  SHELL  version, 
and  the  more  recent  HULL  version  which  treats  the  flow  second  order  accurate 
in  a  Lagrangian  phase  and  then  fluxes  mass,  momentum  and  energy  to  maintain  the 
original  spatial  grid.  The  latter  method  is  the  one  described  in  the  sequel. 
Additional  options  available  include  seven  different  rezone  techniques,  three 
equations  of  state,  single  or  multi -materials  (any  number),  variable  zone  size, 
three  atmospheres,  massless  trace  particles  or  massive  noninteractive  and/or 
interactive  particles  (e.g.,  dust),  water  vapor,  rigid  bodies  within  the  mesh, 
equilibrium  radiation  diffusion,  a  strength  formulation,  magnetohydrodynamics 
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formulation,  and  an  option  for  hand . Ir.g  a  nuclear  air  blast  precursor.  The 
magnetohydrodynamics  formulation  has  recently  been  made  operational  and  the 
strength  formulation  is  still  being  checked  out. 

The  HULL  code  is  based  upon  numerical  techniques  and  algorithms  that  have 
been  developed  by  many  individuals  and  organizations  over  a  period  of  many 
years. 

Historically,  the  development  of  large  Eulerian  hydrocodes  began  with 
Harlow's  SHELL-PIC  code  at  Los  Alamos  in  1955  (ref.  1).  SHELL-PIC  used  the 
particle-in-cell  method  for  the  numerical  solution  of  the  hydrodynamic  equa¬ 
tions.  Due  to  its  Eulerian  nature,  SHELL-PIC  was  able  to  effect  a  solution  to 
the  hydrodynamic  equations  even  when  large  distortions  were  present  in  the 
fluid. 

To  provide  a  good  theoretical  picture  of  the  phenomenology  of  atmospheric 
detonations,  the  AFWL  acquired  two  versions  of  SHELL,  developed  at  General 
Atomic  by  Freeman,  Johnson  and  others  for  hypervelocity  impact  studies.  The 
first  version,  SHELL-PIC,  was  based  on  Harlow's  code.  The  second  code, 

SHELL-OIL,  was  a  continuum  version  (ref.  2).  Comparison  of  PIC  and  OIL  calcula¬ 
tions  of  late  time  phenomenology  of  the  BLUEGILL  Event  (ref.  3)  demonstrated 
the  superiority  of  the  continuous,  SHELL-OIL  method  over  the  discrete,  particle- 
in-cell  method  for  that  class  of  problem. 

In  1965,  at  the  request  of  Defense  Atomic  Support  Agency,  the  AFWL  modified 
SHELL-OIL  (hereafter  called  SHELL)  to  allow  the  calculation  of  high-explosive 
detonation  phenomena  (ref.  4).  The  first  event  calculated  was  DISTANT  PLAIN 
(refs.  5  and  6).  The  AFWL  then  added  a  multi -material  capability  to  SHELL  to 
allow  preshot  computation  of  DISTANT  PLAIN  Event  6  (ref.  7).  Since  then,  the 
AFWL  has  continued  to  provide  numerical  calculations  in  support  of  Defense 
Nuclear  Agency's  high  explosive  tests  (refs.  8  through  12). 

In  addition  to  the  high  explosive  computations,  the  AFWL  performed  over  40 
calculations  of  various  yields  and  heights  of  burst  during  the  period  1965  to 
1970  in  support  of  systems  vulnerability  and  ionization  threat  analyses  (refs. 
13  and  14).  During  this  same  period,  SHELL  was  modified  to  study  dust  and 
pebble  clouds  (refs.  15,  16  and  17).  The  incorporation  of  variable  zone  size 
permitted  the  study  of  nuclear  air  blast  precursors  (refs.  18  and  19).  The 
addition  of  fine  material  capability  and  radiation  diffusion  enabled  us  to 
study  nuclear  weapon  detonation  phenomena  (refs.  20  and  21). 
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Since  1970,  other  versions  of  SHELL  were  written  in  spherical  (ref.  22) 
and  Cartesian  (ref.  23)  coordinates.  SHELL  was  also  rewritten  to  run  on  the 
ILL IAC  IV  computer  (ref.  24). 

By  1971  many  versions  of  SHELL  existed  (more  than  20  separate  decks)  making 
modifications  to  the  code  difficult.  This  situation,  when  coupled  with  the 
desire  to  improve  the  finite  difference  scheme  in  SHELL,  led  Durrett  and  Matuska 
of  AFWL  to  develop  the  HULL  system.  In  order  to  contain  in  one  system  a  great 
multiplicity  of  versions,  it  was  decided  to  write  HULL  in  a  higher  ordered 
language.  This  language  required  a  compiler,  which  was  written  by  Captain  Lewis 
Gaby  of  AFWL  and  was  given  the  name  SAIL.  The  output  of  SAIL  is  primarily 
FORTRAN,  although  some  assembly  code  does  exist.  With  the  higher  ordered 
language,  Durrett  and  Matuska  were  able  to  structure  HULL  so  that  (1)  modifica¬ 
tions  to  the  various  versions  of  the  code  could  be  made  easily,  (2)  typical 
calculations  were  run  three  times  as  fast,  and  (3)  restrictions  on  the  number 
of  zones  or  number  of  materials,  etc.,  disappeared  (i.e.,  a  problem's  size 
becomes  limited  only  by  the  computer's  capabilities  in  terms  of  central  memory, 
disc  space,  etc.,  instead  of  DIMENSION  statements). 

Since  1971  additional  modifications  have  been  made  in  the  HULL  system  by 
other  contributors  such  as  Clemens  of  SAI  (advanced  plotting  and  system  capa¬ 
bilities),  Fry  of  AFWL  (magnetohydrodynamics).  Tower  of  AFWL  (strength  formula¬ 
tion),  Gaby  of  AFWL  (radiation  diffusion),  Chambers  of  AFWL  (nuclear  air  blast 
precursors),  Needham  of  AFWL  (advances  in  burn  routines  and  equations  of  state), 
Prentice  of  AFWL  (auxiliary  debugging  aids),  and  Westmoreland  of  AFWL  (extensive 
revisions  to  KEEL  and  one-dimensional  versions  of  HULL). 

By  June  1976,  HULL  was  on  version  88  and  development  continues  for  various 
types  of  problems  as  the  need  arises.  Currently  a  library  of  almost  300  prob¬ 
lems  exists. 

It  is  worthwhile  to  note  that  the  development  for  HULL  has  primarily  been 
based  on  the  study  of  fireball  growth  and  rise  phenomena,  and  the  prediction 
of  free  field  air  blast  as  well  as  the  effect  of  rigid  structures  on  the  air 
blast.  For  these  applications,  HULL  has  been  shown  to  perform  very  well  when 
utilized  in  the  fashion  that  was  intended. 

The  remainder  of  this  report  discusses  the  architecture  and  operating 
philosophy,  the  differential  equations  and  the  difference  technique,  and 
experimental  and  theoretical  comparisons. 
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SECTION  II 

ARCHITECTURE  AND  OPERATING  PHILOSOPHY 

This  section  briefly  describes  the  architecture  and  operating  philosophy  of 
HULL.  The  motivation  for  developing  HULL  was  to  provide  (1)  flexibility  in 
making  revisions  either  to  improve  or  add  options,  (2)  speed  while  executing 
and  (3)  removal  of  restrictions  on  size  that  are  dictated  by  the  inadequacies  of 
pure  FORTRAN.  In  order  to  accomplish  these  goals  it  was  necessary  to  write 
HULL  in  a  language  of  higher  level  than  FORTRAN;  this  is  referred  to  as  the 
Executive  Language. 

It  was  also  necessary  to  write  a  compiler  for  the  Executive  Language,  which 
was  named  SAIL.  The  combination  of  Executive  Language,  F0RTRm?«,  and  assembly 
language  that  comprises  the  HULL  system  is  a  SAIL  input  file.  It  consists  of 
a  set  of  card-images  which  the  user  can  modify  if  and  when  appropriate.  With 
the  flexibility  obtained  with  the  higher-level  language  it  was  possible  to  auto¬ 
mate  many  burden  tasks,  but  at  the  expense  of  many  man-years  of  code  development 
closely  tied  to  the  operating  system  being  utilized.  With  this  generality 
available  in  the  HULL  system,  auxiliary  programs  were  developed  to  perform 
some  specialized  functions,  many  of  which  were  designed  to  automate  as  much 
as  practical  the  use  of  HULL  on  a  day  to  day  basis.  These  include  programs 
BOW,  PLANK  and  STERN. 

BOW  creates  and  maintains  a  library  tape  (or  LIBTAPE)  which  contains  infor¬ 
mation  on  the  status  of  various  jobs  or  problems  running  under  the  HULL  system. 
Each  problem  run  by  HULL  is  assigned  a  unique  problem  number  for  identification. 
BOW  basically  keeps  a  current  record  of  the  data  tapes  for  each  problem,  and 
has  a  list  of  unused  tapes  which  are  assigned  to  the  various  problems  as 
necessary.  STERN  is  executed  at  the  end  of  each  HULL  run  and  updates  the 
library  tape  with  information  on  the  tapes  used  during  the  run. 

PLANK  builds  the  HULL  code  for  each  run  using  as  input  cards  read  from  the 
input  deck,  information  from  the  HULL  data  tape,  and  the  card  image  information 
on  the  launch  file. 

To  gain  insight  into  how  these  are  utilized,  the  following  example  shows 
how  a  typical  run  might  proceed. 
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The  INPUT  file  is  read  into  the  computer.  It  consists  of  control  cards  and 
other  input,  generally  in  the  form  of  a  card  deck.  Execution  of  this  run  begins 
by  fetching  the  HULL  code  from  the  disc.  Then,  user-generated  changes  are  made 
to  the  codes  that  are  specified  in  the  INPUT  file,  and  the  remainder  of  the 
INPUT  file  is  copied  to  TAPE5.  The  library  tape  is  then  assigned  to  the  job 
and  BOW  executes.  The  BOW  code  reads  the  problem  numbers  and  other  start  param¬ 
eters  from  TAPE5.  It  then  reads  the  library  tape,  finds  the  appropriate  problem 
record,  and  selects  the  correct  HULL  data  tape  to  restart  the  calculation.  BOW 
then  asks  the  operator  to  mount  the  HULL  data  tape,  which  contains  all  the 
hydrodynamic  variables,  by  reference  to  its  location  in  the  computer  tape 
library.  When  this  tape  is  mounted,  BOW  terminates. 

PLANK  executes  next.  It  reads  the  problem  number  and  other  start  informa¬ 
tion  from  TAPE5.  It  then  reads  the  HULL  data  tape  and  selects  the  appropriate 
cycle  for  restart.  Having  found  the  correct  cycle,  PLANK  reads  certain  infor¬ 
mation  off  the  tape  and  with  this  information  decides  which  HULL  code  to 
generate.  PLANK  is,  in  addition  to  other  things,  a  compiler  which  reads  the 
Executive  Language  from  the  card  images  on  the  SAIL  file  and  creates  from  them 
a  complete  FORTRAN  code  which  it  puts  cn  an  appropriate  file.  At  this  point 
PLANK  terminates.  The  FTN  compiler  is  then  invoked  and  compiles  the  FORTRAN 
code.  The  compiled  version  of  HULL  is  then  loaded  and  executed,  and  will  run 
until  a  stoo  condition  is  sensed.  This  is  typically  4  to  5  hours  or  until 
problem  completion. 

STERN  executes  after  HULL  terminates.  STERN  requests  that  the  library  tape 
be  mounted  and  then  updates  the  problem  information  record  to  reflect  the 
changes  made  in  problem  status  on  this  run. 

There  also  exists  in  the  HULL  system  great  amounts  of  code  that  make  it 
possible  to  safely  allow  operators  to  terminate  calculations,  and  then  to  later 
restart  these  calculations  at  the  point  where  they  were  previously  terminated, 
without  any  user  intervention.  The  operator  only  need  resubmit  the  same  deck. 
This  automatic  restart  feature  has  evolved  at  the  AFWL  during  the  last  two 
decades.  This  feature  makes  life  somewhat  more  tolerable  during  periods  of 
frequent  system  crashes,  power  outages  and  other  abnormal  means  of  program 
termination. 

The  data  tapes  that  are  written  and  contain  the  data  necessary  to  restart 
serve  as  a  large  data  bank  for  future  analysis.  Since  each  of  these  calculations 
take  tens  of  hours  on  a  CDC  6600,  it  is  desirable  to  save  ;he  data. 
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SECTION  III 

DIFFERENTIAL  EQUATIONS  AND  DIFFERENCE  TECHNIQUES 

This  section  presents  the  differential  equations  being  solved  with  HULL 
and  the  difference  techniques  utilized.  This  report  considers  only  the  hydro- 
dynamic  calculation.  Other  reports  will  address  the  magnetohydrodynamics, 
radiation  diffusion,  and  the  strength  formulations. 

1.  THE  HULL  2D  DIFFERENCE  METHOD 

The  difference  method  employed  by  the  HULL  code  is  a  fully  second  order 
accurate  method  developed  by  Matuska.  The  equations  below  describe  a  compres¬ 
sible,  nonconducting  and  inviscid  fluid. 

Conservation  of  Mass 

p(V-u)  -  0  (1) 

Conservation  of  Momentum 

p  du  +  Vp  =  -  (2) 

Conservation  of  Energy 

p  +  V* (pu)  =  -  pu*g  (3) 

Equation  of  State 

P  •  f  (p,I)  (4) 


where 

p  =  material  density  (gm/cm3) 

P  =  pressure  (dynes/cm2) 
u  =  (u,v)  fluid  velocity  (cm/sec) 

E  =  I  +  (ergs/gm) 
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I  =  internal  specific  energy  (ergs/p) 
g  =  gravitational  acceleration  (cm/gm2) 
t  =  time  (sec) 

The  procedure  is  to  solve  the  first  four  equations  in  two  phases.  In  phase 
I,  we  solve  the  equations  within  a  Lagrangian  framework.  Equations  (2)  and  (3) 
are  treated  in  this  phase.  The  second  phase  treats  equation  (1)  in  a  manner 
which  causes  dissipation  of  kinetic  energy  into  internal  energy,  especially  in 
regions  of  large  velocity  gradients. 


dV  , 
p  it  + 


=  -pg 


dE  .  1  a(RPU)  .  3  fpvx  , 
p3t  +  R  3R +  9 T(PV) 


where 

R  =  radial  coordinate 
Z  =  axial  coordinate 

U  =  component  of  u  in  the  radial  direction 
V  =  component  of  u  in  the  axial  direction 

In  establishing  finite  difference  analogs  to  equations  5  through  8  we  con¬ 
sider  a  discrete  subset  of  F(R,Z,T)  by  defining 

F(I,J,N)  =  F(R(.I),Z(J),T(N)) 

where  R(I),  Z(J),  and  T(N)  are  particular  values  of  R,  Z,  and  T,  respectively, 
and  the  I,  J,  and  N  assume  integer  values  in  the  range  1  to  IMAX  for  I,  1  to 
JMAX  for  J,  and  0  to  NMAX  for  N.  The  R(I)  and  I ( J )  are  defined  in  terms  of  a 
given  set  of  DR( I )  and  DZ(J)  such  that 
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R ( I )  =  R ( 0 )  +  (SUM,  K=1 ,  1-1,  (DR( K) ) )  +  DR(I)/Z  for  1=2 . IMAX 

R(l)  =  R (0)  +  DR ( 1 )/2 

R(J)  =  Z(0)  +  (SUM ,K=1 ,J-1 , (DZ(K) ) )+DZ( J)/2  for  J-2 . JMAX 

Z(l)  =  Z(0)  +  DZ(l)/2 

where  R(0)  and  Z( 0)  have  some  specified  values. 

The  hydrodynamic  variables  p,  U,  V,  and  I  (internal  specific  energy)  are 
defined  for  each  set  of  coordinates  (I,J)  at  a  particular  time  T(N) .  The 
pressure  P(I,J,N)  is  defined  at  each  point  by  the  equation  of  state  (equa¬ 
tion  4). 

Interpolated  values  for  the  hydrodynamic  variables  of  the  form 
F(I+1/2,J,N),  F ( I ,J+1/2,N) ,  or  F(I,J,N+l/2),  or  similar  forms,  are  de¬ 
fined  in  terms  of  the  F(I,J,N).  In  general, 


F( 1+1/2 ,J ,N)  =  (F(I+1,J,N)  +  F(I  ,J,N) )/2 


F(I,J+1/2,N)  =  (F(I,J+1,N)  +  F(I,J,N))/2. 

This  definition  will  apply  except  where  explicitly  noted, 
a.  Phase  I 

In  method  2  the  finite  difference  analogs  to  equations  (6)  through 
(8)  are  chosen  as 

U ( I , J ,N+1  )  =  U(I,J,N)-DT*(P(I+l/2,J,N+l/2)-P(I-l/2,J,N+l/2)) 

/(p(I,J,N)*DR(I))  (9) 


V ( I , J ,N+1 )  = 


E(I,J,N+1 )  = 


V(I,J,N)-DT*(P(I ,J+l/2 ,N+l/2)-P(I,J-l/2,J,N+l/2)) 
/(p(I,J,N)*DZ(J) )-DT*G( J)  (10) 

E(I,J,N)-DT/p(I ,J,N)*( (R(I+l/2)*P(I+l/2,J,N+l/2) 
*U(I+l/2,J,N+l/2)-R( 1-1/2 )*P( 1-1/2, J.N+1/2) 
*U(I-l/2,J,N+l/2))/(R(I)*DR(I))+(P(I,J+l/2,N+I/2) 
*V(I,J+l/2,N+l/2)-P(I ,J-l/2,N+l/2)*V(I,J-l/2,N+l/2) ) 
/DZ(J)/-DT*V( I ,J ,N+1 )*G(J)  (11) 


'.v.  ' a, 


where 

DT  =  T(N+1)-T(N) 

R( 1+1/2)  =  R(I)+DR(I)/2 
R( 1-1/2)  =  R ( I ) -DR( I )/2 
Z(J+l/2)  =  Z(J)+DZ(J)/2 
Z(J-l/2)  =  Z(J)-DZ(J)/2 
G(J)  =  VALUE  OF  G  AT  Z(J). 

All  the  values  appearing  in  equations  (9)  through  (11)  are  imnediately 
known  except  the  time  advanced  (N+l/2)  values  for  pressure  and  velocity.  These 
time  advanced  values  are  used  so  that  the  approximations  to  the  partial  deriva¬ 
tives  appearing  in  equations  (6)  through  (8)  may  be  centered  in  time  and  space. 
In  the  case  where 


and 


DR( I )  =  constant  For  I=1,...,IMAX 

DZ(J)  =  constant  For  J=1 ,. . . ,JMAX, 


this  produces  a  fully  second  order  accurate  difference  method.  In  a  region 
where  the  DR(I)  and  DZ(J)  are  not  constant  the  second  order  accuracy  is  lost. 
This  adversely  affects  the  stability  of  the  first  phase  calculation.  The 
amount  of  instability  which  may  be  obtained  is  related  to  the  magnitude  of  the 
incremental  changes  in  DR(I)  and  DZ(J). 

Most  of  the  computations  in  the  first  phase  are  expended  in  obtaining 
the  time  advanced  values  for  pressure  and  velocity.  The  time  advanced  veloci¬ 
ties  are  obtained  by  differencing  equations  (6)  and  (7)  as 


U( 1+1/2, J, N+l/2)  =  U(I+l/2,J,N)-DT/(2*p(I+l/2,J, N+l/2)) 

*( (P(I+1 ,J,N)-P(I,J»N))/(R( 1+1 )-R( I) ) )  (12) 

V ( I , J+l/2. N+l/2)  =  V(I,J+l/2,N)-DT/(2*p(I,J+l/2, N+l/2)) 

*( (P( I,J+1 ,N)-P(I ,J,N) )/ (Z( J+l ) — Z( J ) ) ) 
-G(J+l/2)*DT/2  03) 


where 


G(J+l/2)  =  ( G ( J ) +G (J+l ) )/2. 
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The  time  advanced  densities  appearing  in  equations  (12)  and  (13)  are  obtained 
by  differencing  equation  (5)  as 

p(I+l/2,J,N+l/2)  =  p(I+l/2,J,N)*(l-DT/(2*R(I+l/2))*(R(I+l) 

*U(I+1 , J ,N)-R( I )*U( I , J,N) )/ (R( 1+1  )-R( I ) ) )  (14) 

p(I,J+l/2,N+l/2)  =  p(I,J+l/2,N)*(l-DT/2*(V(I,J+l,N)-V(I,J,N)) 

/ (Z(J+1  )-Z(J) ) )  (15) 


where 


p(I+l/2,J,N)  =  M(I,J,N)+M(I+L,J,N))/  (tt*(R(  1+3/2  )**2 
-R(I+1/2)**2)*DZ(J) ) 

p(  I , J+1/2.N)  =  (M(I,J,N)+M(I,J+1,N))/tt*(R(I+1/2)**2 
-R(I-1/2)**2)*(DZ'(J)+DZ(J+1 ))) 

and  the  mass  associated  with  a  point  (I.J.N)  is  defined  by 

M(  I ,  J  ,N)  =  p(I,J,N)*(tt*(R(I+1/2)**2-R(I-1/2)**2)*DZ(J)) 

where  tt  =  3.14159...  and  R(  1+3/2)  =  R(  1+1 )  +  DR(I+l)/2. 

The  time  advanced  pressure  appearing  in  equations  (9)  through  (11) 
requires  a  little  more  effort.  First,  an  alternative  energy  equation  can  be 
obtained  from  equations  (2)  and  (3)  as 

p  $  +  •  u  =  0  (16) 

An  effective  y  can  be  defined  by 

T',+5T  (17) 

We  will  assume  for  the  purposes  of  calculating  a  half  time  step  advanced  pres¬ 
sure,  which  in  turn  is  used  in  approximating  the  partial  derivatives  in  equa¬ 
tions  (9)  through  (11),  that  the  Lagrangian  derivative  with  respect  to  time  of 
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Y  is  small  and  can  be  ignored.  In  application,  it  is  only  required  that  the 
change  in  y  at  a  particular  point  be  small  over  a  time  of  DT/2. 

Taking  the  Lagrangian  derivative  with  respect  to  time  in  equation  (17) 
and  using  equations  (1)  and  (16),  we  can  write 


f+YP?,J:0 


(18) 


Equation  (18)  is  used  to  obtain  time  advanced  pressures  given  by 

P(I+1/2,J ,N+1/2)  =  P( 1+1/2 ,J,N)*(l-DT*y( 1+1/2,  J,N)*(R(I+1) 

*U(I+1 ,J,N)-R(I)*U(I,J,N))/(2*R(I+1/2)*(R(I+1) 

-R(I)))) 

P(I,J+l/2,N+l/2)  =  P(I,J+1/2,N)*(1-DT*y(I,J+1/2,N) 

*(V(I,J+1 ,N)-V (I ,J  ,N) )/ (2*(Z(J+l )-Z(J) ) ) ) 

where  y  is  defined  from  equation  (17)  as 

y(I+l/2,J,N)  =  1+P(I+1/2,J ,N)/(p( 1+1/2,  J  ,N)*I (1+1/2, J ,N) ) 

y(I,J+l/2,N)  =  l+P(I,J+l/2 ,N)/(p(I ,J+1/2,N)*I(I ,J+1/2,N) ) . 

All  quantities  needed  to  solve  equations  (9),  (10),  and  (11)  are  now  defined. 
Solution  of  these  equations  will  complete  a  second  order  accurate  Lagrangian 
calculation.  The  next  step  would  normally  be  that  of  transporting  mesh  vertices. 
Instead  we  choose  to  flux  the  hydrodynamic  quantities  to  retain  the  original 
mesh  configuration.  This  calculation  is  done  in  phase  II. 

b.  Phase  II 

Changes  in  density  are  computed  in  phase  II  by  calculating  a  mass  flux 
between  mesh  points  and  then  transporting  the  appropriate  amount  of  mass  from 
point  to  point.  The  transported  mass  carries  with  it  a  proportionate  amount 
of  internal  energy  and  momentum.  The  velocities  and  specific  internal  energy 
are  then  redefined  at  each  mesh  point  by  conserving  momentum  and  total  energy 
at  that  point. 

15 
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The  mass  flux  between  mesh  points  is  defined  as  the  product  of  the 
interpolated  velocity,  the  density  as  defined  by  solution  of  equation  (1), 
the  intermediate  cross  sectional  area,  and  the  time  step.  The  equations  are: 

MF(I+1/2,J,N+1 )  =  U(I+l/2,J,N+L)*p( 1+1/2  ,J,N+1  )*2*tt*R(  1+1/2) 

*DZ(J)*DT 

MF ( I  ,J+1/2,N+1 )  =  V(I,J+1/2,N+1  )*p( I, J+1/2.N+1  )*2*tt*R( I )*DR ( I ) 

*DT 

where  the  time  advanced  densities  are  obtained  by  differencing  equation  (5)  as 

p(I+l/2,J,N+l )  =  p ( ID,  J,N)*(1-DT/R( 1+1/2)* (R( 1+1 ) 

*U(I+1 ,J,N+1 )-R(I)*U(I ,J,N+1 )  )/(R(I+l )-R(I) )) 

p(I,J+l/2,N+l)  -  p(I,JD,N)*(l-DT*(V(I,J+l,N+l)-V(I,J,N+l)) 
/(Z(J+1)-Z(J))1 


where 

ID  =  1  If  U(I+1/2,J,N+1)  GT  0 

=1+1  If  U(I+1/2,J,N+1)  LT  0 

JD  =  J  If  V(I ,J+1/2,N+1 )  GT  0 

■  J+l  If  V(I ,J+1/2,N+1 )  LT  0. 

This  is  the  classical  donor  cell  differencing  technique.  The  most  obvious 
advantages  of  this  technique  are  its  rigid  numerical  conservation  and  its 
stability.  This  scheme  also  insures  that  more  material  cannot  be  removed 
from  a  point  than  is  present. 

HULL  has  a  continuous  rezone  capability.  When  this  is  employed,  the 
interpolated  velocities  appearing  in  the  mass  flux  equations  are  replaced  by 

U(I+1/2,J,N+1)  -  UR(I+1/2,J,N+1) 

V(I,J+1/2,N+1)  -  VR(I,J+1/2,N+1) 
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where  UR  and  VR  are  the  interpolated  grid  velocities  (determined  arbitrarily 
by  how  fast  one  wishes  to  transport  the  coordinate  grid).  The  corresponding 
momentum  fluxes  are 

UF(I+1/2,J,N+1 )  =  MF( 1+1/2 ,J ,N+1 )*U ( ID,J  ,N+1 ) 

VF(I+1/2,J,N+1)  =  MF(I+l/2, J ,N+1 )*V ( ID, J, N+l ) 

UF ( I ,J+1/2,N+1 )  =  MF ( I , J+l/2, N+l )*U  ( I  ,JD,N+1 ) 

VF ( I ,J+1/2,N+1 )  =  MF(I ,J+1/2,N+1 )*V(I,JD,N+1 ) 

and  the  energy  fluxes  are 

EF(I+1/2,J,N+1)  =  MF(I+1/2,J ,N+1 )*E(ID,J,N+1 ) 

EF ( I ,J+1/2,N+1 )  =  MF ( I , J+l/2 ,N+1 )*E( I ,  JD.N+l ) 

When  these  quantities  are  fluxed,  final  values  for  mass,  density, 
velocity,  and  energy  are  computed  by 

M  ( I ,  J )  =  M(I,J,N)+MF(I-1/2,J,N+1)+MF(I,J-1/2,N+1) 
-MF(I+1/2,J,N+1 )-FM(I , J+l/2, N+l) 

p(I,J)  -  M( I ,J)/(tt*(R(1+1/2)**2-R(I-1/2**2)*DZ(J) ) 

U(I,J)  =  (U ( I , J , N+l  )*M(I ,J,N)+UF( 1-1/2 ,J ,N+1 )  +U F ( I , J - 1 /2,N+1 ) 

-UF( I+1/2.J, N+l )-LF( I ,J+l/2 ,N+1 ) )/M( I , J) 

V(I,J)  =  (V ( I , J ,N+1  )*M( I , J ,N )+EF (I-1/2.J, N+l  )+EF(I  ,J-1/2,N+1 ) 

-EF( 1+1/2,  J  ,N+1)-VF(I+1/2,J , N+l ))/M(I,J) 

1(1, J)  =  (E(I , J ,N+1 )*M(I,J,N)+EF(I-1/2,J ,N+1  )+EF(I,J-l/2,N+l ) 

“EF (I+1/2.J ,N+1 )-EF (I.J+1/2 ,N+1 )-(U(I,J)**2+7(I,J)**2) 
*M(I,J)/2)/M(I,J) 

E ( I ,  J )  =  I(I,J)+(U(I,J)**2+V(I,J)**2)/2 
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where  the  lack  of  a  time  specification  indicates  final  values  for  this  time 
step 

2.  HULL:  THE  3-DIMENSIONAL  VERSION 

The  three-dimensional  version  of  the  HULL  code  solves  the  finite  difference 
analogs  of  the  system  of  partial  differential  equations  describing  inviscld, 
nonconducting  fluid  flow  in  the  form 


+  p7  ■  u  ■  0 

(18) 

P  f  r  VP  =  P? 

(19) 

p  +  V  •  Pu"  =  pil  •  g 

(20) 

P  =  (Y-l)pI 

(21) 

where 

p  =  material  density  (gm/cm3) 

P  =  pressure  (dynes/cm2) 
u  =  (u j v jW)  the  fluid  velocity  (cm/sec) 

E  =  I  +  ^  (ergs/gm) 

I  =  internal  specific  energy  (ergs/gm) 
g"  =  acceleration  due  to  gravity  (cm/sec2) 
t  =  time  (sec.) 

Y  =  specific  heat  at  constant  pressure  divided  by  specific  heat 
at  constant  volume 

Equations  (18)  through  (21)  are  solved  in  two  phases.  The  first  phase  pro¬ 
ceeds  as  a  Lagrangian  calculation  insuring  exact  compliance  with  equation  (18) 
when  all  mass  points  are  properly  transported.  Choosing  the  following  notation 
for  the  finite  difference  analogs  of  the  above  equation, 

f(xi’yj’zk,tn)  “  fi,j,k 
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and  where  subscripts  on  the  left  and  right  sides  of  the  equations  are  suppressed 
on  the  right-hand  side,  the  finite  difference  analogs  of  equations  (19)  and  (20) 
in  Cartesian  three  space  with  the  z  axis  oriented  in  the  direction  of  g-  are: 


u?+!  -  un  -  JL  (p£j?  -  pf*) 

pnDX  1  2  12 


i  ,j,k 


(22) 


v?+ !  =  vn  -  -P-  (p"S  -  pl+i) 


i  »j  *k 


pnDY  J-V 


W?lk  -  W"  -  ±  <5  -  O  +  DT'Gk 


p  DZ 


(23) 


E?lk  - E"  -  t  O/0x 


-  <pjn:s  »j^)/w 


+  (C  <5  -  C 


+  DT-V-G, 


(24) 


where 


DT  =  Tn+1  -  Tn 


OX  -  x.+1 


DY 

DZ 


Vl  '  Yj 

zk+l  ■  h 


The  time  and  space  centered  quantities  are  defined  from  the  previous  time 
step  mesh  quantities.  In  addition,  we  need  pressure  and  velocity  terms  in  the 
above  equations  and  use  an  alternate  form  of  the  energy  equation  obtained  from 
equations  (19)  and  (20)  for  these; 


dl 


Ht 


+  PV  •  u  =  0 


(25) 


Solving  equations  (18),  (21),  and  (25)  for  a  constant  gamma,  we  obtain 
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dP 

dt 


TT  +  Y  PV  •  U  =  0 


(26) 


The  assumption  of  a  constant  ganma-law  gas  is  not  critical  to  the  methodology 
since  an  effective  ygff  =  +  1  can  almost  always  be  used  since  y  does  not 

change  significantly  in  a  time  step. 

The  finite  difference  analogs  of  equations  (19)  and  (26)  for  the  time  and 
space  centering  quantities  used  in  equations  (22),  (23),  and  (24)  are: 


Mn+}S 

Ui+Js, j  »k 


+  Ci 


DT 


DX(o'j+J5  +  p 


n+^s 

i+1 


(P" 


+1 


P") 


(27) 


V 


n+h 

i  ,j+}5,k 


*1  - 


DT 


dy(c"+h  + 


_  Ip” 

,n+^\  j+1 

3j+l; 


Wi,j,k+J5 


•  ill  .  i  ■  n 
.  uk  *  wk+i 


DT 


DZ(p"+i5 


+  O 


(P 


k+1 


-  pjj)  +  k+ij 


(28) 


(29) 


pn+h 
i+Js,j  ,k 


Pn+J2 

i  .j+J2»k 


pn+Js 
i  ,j,k+% 


+  pn 

Ki+1 
1 - 


+  P 


itl 


+  p 


k+l 


"  Y  Ox  ^Ui+1 
'  Y  OX  (Vj+l 
\}  "  Y  OZ  ^Wk+1 


(30) 

(31) 

(32) 


The  G  terms  in  equations  (23),  (24),  and  (29)  are  calculated  by  solving  the 
hydrostatic  equations  on  the  finite  difference  mesh  at  the  appropriate  mesh 
interval.  (An  automatic  program  abort  occurs  if  this  calculated  value  of  G 
ever  differs  from  the  correct  theoretical  value  at  the  current  altitude  by 
more  than  1  percent.) 

The  above  set  of  £inite  difference  equations  is  second  order  accurate  for 
a  fixed  mesh  interval  ar.d  is  reasonably  well-behaved  when  adjacent  cells  are 
within  10  percent  of  one  another  in  size. 
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Normally,  a  Lagrangian  solution  of  equations  (18),  (19),  and  (20)  requires 
the  transport  of  vertices.  In  the  three-dimensional  mesh  we  choose  instead  to 
flux  across  cell  boundaries  while  conserving  energy,  momentum  and  mass.  Al¬ 
though  this  is  the  main  source  of  entropy  production  in  the  method,  it  has  the 
advantage  of  providing  sufficient  dissipation  to  permit  calculation  of  strong 
shocks  without  the  difficulty  of  managing  an  explicit  viscosity  term. 

In  addition,  a  continuous  rezone  is  overlayed  on  the  flux  scheme.  In 
three- dimensions  it  is  necessary  to  alternate  the  application  of  rezone 
velocities  so  that  only  one  cell  boundary  is  moved  during  a  given  time  step. 
This  requirement  arises  due  to  accumulation  of  second  and  higher  order  errors 
in  ambient  zones  which  can  upset  hydrostatic  equilibrium. 

By  solving  equation  (18)  we  obtain  the  material  density  to  be  fluxed. 

Taking  account  of  dilatational  terms  in  all  three  directions  simultaneously 
in  a  straightforward  manner  leads  to  instability  in  regions  of  rapid  expansion. 
A  method  of  first  order  accuracy  but  of  superior  stability  characteristics  can 
be  obtained  by  solving  for  the  dilatational  terms  one  coordinate  at  a  time,  and 
since,  if  stability  conditions  are  adhered  to,  the  eigenvalues  of  each  of  the 
finite  difference  operators  are  less  than  unity,  then  their  products  will  also 
be  less  than  one.  This  leads  to 


i+%,j,k  =  pid^Awi 

(33) 

i  >j+}s»k  "  Pjd/Auj 

(34) 

i  ,j  ,k+*s  =  Pkd/Auk 

(35) 

The  d  following  a  cell  index  indicates  tfje  donor  cell  values  are  used.  Also, 


Au, 


i.j.k 


i  +  21  nin+i  nM\ 

1  +Ux  (ui+i  '  V 


i»i  +  21  (vn+1 
1  DY  ^Vj+1 


-  v">]  '  [  +  m' 


(Wn+1V) 
'wk+l  V 

(36) 


The  effect  of  this  apparently  ad  hoc  treatment  is  to  drastically  reduce  the 
rather  bothersome  corner-coupling  difficulties  which  often  result  in  spherical 
features  such  as  fireballs  being  "squared-off."  We  also  found  that  we  can  run 
the  code  at  about  0.9  of  the  usual  one-dimensional  Courant  condition,  an 
obvious  computational  speed  advantage. 
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SECTION  IV 

EXPERIMENTAL  AND  THEORETICAL  COMPARISONS 

This  section  compares  some  theoretical  calculations  using  HULL  (and  its 
predecessor,  SHELL)  with  various  high  explosive  (HE)  experimental  data.  The 
AFWL  has  provided  air  blast  calculations  in  support  of  nearly  all  large  scale 
HE  detonations  beginning  with  the  DISTANT  PLAIN  series  and  continuing  through 
dice  THROW  (refs.  4  through  14).  Hence,  there  exists  a  wealth  of  data  for  com- 
pari  sons  with  theoretical  calculations. 

Figure  1  is  a  comparison  of  the  overpressure  versus  ground  range  for  a 
500  ton  tangent  sphere  of  TNT.  This  calculation  was  made  with  our  vision  o 
the  SHELL  code  (ref.  6),  which  now  has  been  replaced  by  the  HULL  code  at  the 
WL  Note  that  although  SHELL  results  fall  below  the  experimental  data  below 
pressures  of  about  10  psi ,  HULL  results  are  considerably  improved  over  those 
from  SHELL  in  that  pressure  regime.  Experimental  results  are  plotted  rom 
PRARIE  FLAT  (ref.  9)  and  DIAL  PACK  which  were  both  detonated  at  Suffiel  , 

Canada. 

MINE  UNDER  (refs.  5  and  10)  was  a  100  ton  TNT  sphere  detonated  1  diameter 
above  the  ground.  Agreement  for  the  peak  overpressure  versus  ground  range 
between  calculation  and  experiment  was  similar  to  that  of  the  ta-.gent  sp  ere 
case.  Positive  phase  duration  (fig.  2)  is  somewhat  more  difficult  to  measure 
experimentally  than  is  peak  overpressure.  This  is  reflected  in  the  relatively 

large  experimental  error  bars. 

Figures  3  through  6  compare  several  parameters  measured  during  MIXED  COMPANY, 
a  500  ton  tangent  sphere  of  TNT  detonated  near  Grand  Junction,  don* .  with 
AFWL  calculations.  Measurements  were  made  by  Ballistics  Research  Laboratory 
(Brl)  and  the  Civil  Engineering  Division  of  the  AFWL.  Calculations  were  made 

by  the  AFWL  Technology  Division. 

The  peak  overpressure  is  in  good  agreement  as  usual.  The  overpressure 
impulse  (fig.  4),  although  somewhat  more  difficult  to  measure,  appears  to  be 
in  good  agreement.  As  mentioned  previously  the  positive  duration  is  very 
difficult  to  measure,  and  is  reflected  by  the  large  amount  of  experimental 

scatter. 
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COMPANY  (positive  phase  duration) 


ConsideraDle  comment  has  been  made  regarding  the  theoretical  capability  to 
properly  calculate  dynamic  pressure.  From  the  conservation  laws  of  mass,  momen¬ 
tum  and  energy  one  can  infer  that  if  one  parameter  at  the  shock  front  is  known, 
all  other  parameters  are  uniquely  determined.  Since  the  codes  solve  the  con¬ 
servation  equations,  it  is  argued  that  when  the  overpressure  is  being  accurately 
calculated,  the  dynamic  pressure  must  also  be  accurately  calculated. 

Because  dynamic  pressure  cannot  be  measured  directly  with  present  instrumen¬ 
tation,  the  method  used  is  to  measure  the  total  pressure  and  the  overpressure 
at  the  same  point.  The  dynamic  pressure  is  calculated  by  making  a  point  by 
point  comparison  of  the  two  measured  waveforms.  This  can  lead  to  accentuated 
differences;  and,  therefore,  several  smoothing  techniques  are  used  to  attempt 
to  damp  resulting  oscillations. 

In  addition  to  the  usual  difficulties  mentioned  above,  the  MIXED  COMPANY 
(ref.  12)  total  pressure  gages  were  3  feet  above  ground,  whereas  the  overpres¬ 
sure  was  measured  at  ground  level.  The  calculations  indicate  a  significant 
pressure  gradient  between  ground  level  and  3  feet.  This  would  lead  to  lower 
"measured"  dynamic  pressure.  Figure  6  shows,  in  general,  this  lowering  of  the 
peak  dynamic  pressure. 

Figures  7  through  14  compare  calculated  and  experimental  results  for  Dipole 
West  (ref.  26),  shots  8  through  11.  These  experiments  used  two  charges  deto¬ 
nated  one  above  the  other  such  that  the  distance  between  charges  was  twice 
that  of  the  lower  charge  above  the  ground.  Thus,  direct  comparisons  could  be 
made  between  real  and  ideal  surfaces  and  calculations.  Shots  8  and  11  had  a 
height  of  burst  of  25  feet  (50  feet  charge  separation),  the  ground  surface 
was  smooth  for  shot  8  and  had  approximately  1  foot  corrugations  in  shot  11. 

The  experimental  data  in  figure  7  are  for  the  real  surface  only;  the  calculation 
falls  between  the  arrival  times  for  the  rough  and  smooth  surfaces. 

Figure  8  is  a  comparison  of  the  overpressure  peaks.  The  calculation  falls 
low  by  about  10  percent  at  ground  zero  but  by  a  distance  of  10  feet  has  regained 
the  peak  and  falls  within  experimental  scatter.  These  calculations  were  run 
with  the  HULL  code;  note  that  the  overpressure  peak  agreement  at  10  psi  is 
excellent  when  compared  to  the  smooth  surface.  The  rough  surface  has  slightly 
enhanced  the  overpressure  peak  at  low  overpressure. 

Figure  9  shows  the  effect  of  decreasing  the  height  of  burst  from  25  to  15 
feet.  The  flattening  of  the  curve  near  100  feet  is  caused  by  the  reflected 
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Figure  9. 


Comparison  of  Calculated  and  Experimental  Maximum 
Overpressures  at  Ground  L«vel ;  Shots  9  and  10, 
Dipole  West,  (extracted  from  ref.  26,  p.  213) 
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shock  from  the  second  burst  catching  that  of  the  first.  Again,  the  calculation 
falls  within  5  percent  of  the  experimental  curves. 

Figures  10  and  11  are  for  gages  mounted  3,  10,  and  25  feet  above  the  ground 
surface.  This  demonstrates  that  the  calculations  are  handling  the  entire  flow 
field  properly,  not  just  a  small  region  near  the  ground. 

Figure  12  clearly  shows  the  slowing  of  the  shock  front  by  the  real  surface 
compared  to  the  ideal  used  by  the  calculation.  The  maximum  difference  is  about 
2  percent. 

Dipole  West  (ref.  26)  measurements  for  dynamic  pressure  were  also  made.  The 
significant  difference  here  is  that  the  total  pressure  and  overpressure  were 
measured  at  the  same  ground  range  and  at  the  same  height  above  the  ground.  The 
peak  dynamic  pressure  data  (fig.  13)  has  a  relatively  large  amount  of  scatter 
for  the  reasons  mentioned  earlier.  The  integrated  dynamic  pressure  has  the 
effect  of  smoothing  individual  oscillations,  thus  reducing  the  scatter  and 
agreeing  much  better  with  the  calculations  (fig.  14). 

Some  of  the  most  carefully  performed  experiments  in  recent  years  were  made 
by  Carpenter  at  TRW  (ref.  27).  A  large  number  of  triply  redundant  detonations 
of  PBX-9404  spheres  were  exploded  over  a  polished  concrete  slab. 

Calculations  were  made  by  the  Technology  Division  of  the  AFWL  in  support 
of  this  project.  Figures  15  and  16  show  that  although  the  code  is  missing  the 
absolute  peaks  the  integrated  waveforms  (impulse)  show  excellent  agreement. 

The  most  recent  HE  tests  have  been  those  of  PRE  DICE  THROW.  Event  2  was  a 
6  ton  detonation  of  a  capped  cylinder  of  AN/FO.  The  HULL  calculation  required 
not  only  a  never  before  calculated  geometry  but  a  relatively  unknown  explosive. 
Figure  17  is  a  comparison  of  the  experimental  and  calculated  overpressure  peaks. 
Agreement  is  excellent.  In  this  particular  experiment,  waveform  comparison  of 
actual  measured  data  was  possible.  The  comparison  showed  that  the  calculated 
waveforms  were  virtual  overlays  of  those  measured. 

It  has  been  the  intent  of  this  section  to  acquaint  the  reader  with  the 
accuracy  of  the  methods  employed  in  theoretical  air-blast  calculations  performed 
at  the  AFWL  by  presenting  a  wide  variety  of  experimental  data  compared  with 
calculated  results.  In  nearly  all  cases  discrepancies  between  calculations 
and  experiments  can  be  understood  in  terms  of  minor  shortcomings  of  the 
experimental  or  calculational  technique.  In  no  case  are  any  of  the  remaining 
discrepancies  considered  serious. 
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Figure  17.  PRE  DICE  THROW  1  Calculation  (Event  2) 
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The  large  variety  of  calculations  of  air  blast  indicates  that  the  large 
scale  hydrocode  is  a  reliable  method  of  predicting  all  phases  of  shock  propaga¬ 
tion  from  source  to  a  few  psi.  The  HULL  code  has  been  used  to  calculate  the 
blast  phenomenology  from  high  explosive  sources  ranging  from  8  pounds  to  500 
tons.  With  few  exceptions  the  calculational  results  are  in  excellent  agreement 
with  experimental  data.  Calculations  of  blast  from  nuclear  detonations  are 
somewhat  easier  to  make  than  those  involving  high  explosives  because  the  large 
mass  of  detonation  products  is  not  present.  The  calculation  must  be  concerned 
with  only  one  material,  air,  which  is  well  understood.  On  this  basis,  we 
expect  the  results  of  nuclear  blast  calculations  to  be  at  least  as  accurate 
(if  not  more  so)  as  those  for  high  explosives. 
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Sandia  Lab  POB  969 
Livermore ,  CA  94550 
1  Org  8000,  Dr.  Cook 

1  R.  Scholer 

1  T.  Gold 

1  J.  Mansfield 

LLL,  POB  808  Livermore,  CA  94550 


NRL,  Wash,  DC  20390  1 

Code  2070  1 

J.  Boris  1 

NSWC,  Wht  Oak,  Slvr  Sp,  MD  20910 

Code  730  1 

J.  Peters  2 

Exp  Eff  Div,  Mr.  Hornig 

DNA,  Wash,  DC  20305 

DDST ,  Dr.  Atkins,  Mr  P  Haas  1 

STTL  i  1 


,  S?SS i  '.E.t  Sevin;1  Capt  J.  Stockton; 

Maj  T.  Strong;  K.  Goering 
Dr.  Ulrich 

SPAS,  J.  Moulton;  M.  Rubenstein,  1 
Capt  Garrison,  Maj  Anderson  1 

STAP  1 

RAAE,  P.  Fleming;  Capt  Mueller, 

Maj  Bigoni 

RATN,  Cmdr  Alderson 

8 

ODDR&E/Asst  Dir,  Strat  Wpns, 

Wash  D  C  20301 


L-51,  M.  Hanson 
TID 

Jeff  Thompson 

LASL,  POB  1663,  Los  Alamos  NM  87544 
Rpt  Lib 

J-10,  E.  Jones;  R.  Gentry' 

Aerospace  Corp  POB  92957,  LA,  CA 
90009 
Lib 

R.L.  Beck 

Bm  Eng  Co,  300  Sparkman 
Huntsville,  AL  35807 
V.  Kenner,  MS-44 
C.  Wayne 

J.  Dobkins 

McDon  Doug,  5301  Balso 
Huntington  Beach,  CA  92547 

K.  Stone;  P.  Lewis;  D.  Dean;  T. 
Thomas;  D.  Hildebrand;  J.  Logan; 
R.  Zach,  H.  Herdman 
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No .  cys 


Kaman  Sci,  CO  Spgs,  CO  80907 
D.  Pirio;  D.  Sachs 


Kaman  Av/J.  Reutenik,  83  2nd  Ave 
Burlington,  MA  01803 

MartMar,  POB  5837,  Orlando,  FL  32805 
R.  Heffner;  C.  Washburn; 

H.  Suguichi;  D.  Cooke 

MRC,  1  Presidio,  Sta  Barbara,  CA  93101 
R.  Christian;  D.  Sowle 

SAI,  POB  2351,  La  Jolla,  CA  92037 
C.  Busch;  R.  Lowen 

SAI/D.  Hove,  101  Cont  Bldg 
Ste  310,  El  Segundo,  CA  90245 

SAI/R.  Hillendahl,  POB  10268 
Palo  Alto,  CA  94303 

SAI/M.  Tobrinen,  Ste  808 

1901  N.  Ft  Myer,  Arlington,  VA  22209 

SAI/R.  Clemens,  122  LaVeta  Dr  NE 
Albuquerque  NM  87108 

BTL/W.  Troutman,  Whippany,  NJ  07981 

Calspan/R.  Deliberis,  POB  235 
Buffalo,  NY  14221 

ISI/W.  Dudziak,  123  W.  Padre 
Sta  Barbara,  CA  93105 

Lockheed,  3251  Hanover 
Palo  Alto,  CA  94304 

MITRE,  Box  208  Bedford,  MA  01730 
S.  Morin;  J.  Brown 

CRT/M.  Rosenblatt,  19720  Ventura  Ste  H 
Woodland  Hills,  CA  91364 

R&D  Assoc,  POB  3580,  Sta  Monica  CA  90403 
H.  Brode;  C.  Knowles;  J.  Carpenter 
R.  Lelevier 

MARA/S.  Kahalas,  385  Elliot, 

Newton  Upper  Falls,  MA  02164 

48 


TRW  Sys  Gp,  1  Sp  Pk 
Redondo  Beach,  CA  90278 

1  Lib 

2  A.  Kuhl 

1  J •  Peterson 

1  AVC0/G.  Grant,  201  Lowell 

Wilmington,  MA  01887 

1  IDA,  400  Arm-Nav  Dr 

Arlington  VA  22202 

1  ASD/I&L ,  Wash,  DC  20301 

1  ARM  For  Stf  Col/Lib, 

Norfolk,  VA  23511 

1  .  Ind  Col  Arm  For, 

Ft  McNair,  Wash,  DC  20315 

1  Nat  War  Col/Lib, 

Ft  McNair,  Wash,  DC  20315 

1  USAMC/ AMCRD-RP-B  Wash,  DC  20315 

1  USAEC/AMSEL-RD 

Ft  Monmouth,  NJ  07703 

USNCEL,  Pt  Hueneme,  CA  93041 
1  Code  L31 

1  CEC  Ofcr 

1  MIT/Dr.  Hansen,  77  MA  Ave 

Cambridge,  MA  02139 

1  Princeton  U/PPL,  Dr  Bleakney 

Princeton,  NJ  08540 

Un  MI,  POB  618,  Ann  Arbor,  MI 
48104 

1  Rsch  Soc  Ofc 

1  IS&T,  G.  Frantti 

1  St  Louis  Un/Dr.  Kisslinger 

221  N.  Grand,  St  Louis,  MO  63100 

1  Boeing/R.  Carlson,  Seattle, 

WA  98124 

1  GATC/MRD  Div  Lib,  7501  N  Natchez 

Niles,  IL  60648 
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No.  cys 

1  H&N  Sp  Pro  D±v/S .  Smith 

849  S  Bwdy ,  LA,  CA  90014 

1  URS,  1811  Trousdale 

Burlingame,  CA  94010 

1  SWRI/M.  Whitfield,  8500  Culebra 

San  Antonio,  TX  78228 

1  Un  IL/Dr.  Newmark,  Tal  Lab 

R207,  Urbana,  IL  61803 

1  IITRI,  10  W.  35 

Chicago,  IL  60616 

1  RAC/Lib,  McLean,  VA  22101 

1  ADC  (ADS WO) ,  Kirtland  AFB, 

NM  87117,  Stop  70 

1  .  SACLO ,  Kirtland  AFB,  NM 

87117,  Stop  68 


TACLOS,  Kirtland  AFB,  NM 
87117,  Stop  67 


AFWL,  Kirtland  AFB,  NM  87117 


1 

HO /Dr.  Minge 

2 

SUL 

1 

DE 

1 

DEV 

1 

DEX 

1 

DY 

1 

DYS 

1 

DYV 

10 

DYT 

1 

DYX 

20 

Capt  M.  Stucker/DYT 

20 

Capt  M.  Fry/DYT 

1 

SA 

1 

SAB 

1 

SAS 

2 

SAT/Mr.  Sharp,  Mr.  Mu: 

1 

SAW 

1 

ESD  (DEE),  Hanscom  AFB, 
MA  01730 

1 

Official  Record  Copy 
Capt  Mark  Fry/DYT 
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